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Abstract Statistical analysis of Diffusion Tensor Imaging (DTI) data requires a 
computational framework that is both numerically tractable (to account for the 
high dimensional nature of the data) and geometric (to account for the nonlinear 
nature of diffusion tensors). Building upon earlier studies that have shown that 
a Riemannian framework is appropriate to address these challenges, the present 
paper proposes a novel metric and an accompanying computational framework for 
DTI data processing. The proposed metric retains the geometry and the compu- 
tational tractability of earlier methods grounded in the afhne invariant metric. In 
addition, and in contrast to earlier methods, it provides an interpolation method 
which preserves anisotropy, a central information carried by diffusion tensor data. 

Keywords diffusion tensor MRI • interpolation • spectral decomposition • 
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1 Introduction 



Diffusion-weighted imaging (DWI) allows non-invasive quantification of the self 
diffusion of water in vivo. In biological tissues, characterized by cell membranes 
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and cytostructures, the movement of water is restricted because of these barriers. 
In tissues as white matter, which is highly directional, the resulting movement of 
water is therefore anisotropic. In this way, high diffusion anisotropy reflects the 
underlying highly directional arrangement of white matter flbre bundles. Diffusion 
measurements (which use the same tools as magnetic resonance imaging (MRI)) 
can characterize this anisotropy. The most common representation of this direc- 
tional diffusion is through the us e of diffusion tensors, a formalism introduced by 
Basser et al in 1994 (Basser et all , [l994) . Since then , othe r highe r level representa- 
tions have been int roduced, as the Q-Ball Ima ging ( Tuchl . [2Qo3 ) and the Diffusion 



Kurtosis Imaging ( Jensen and Helpernl . |2Q1Q[ ). In the context of Diffusion Tensor 



Imaging (DTI), each voxel of the image contains a diffusion tensor, which is de- 
rived from a set of DWI measured in different directions. A diffusion tensor is 
nothing else than a symmetric positive deflnite matrix whose general form is given 
by 

(Dxx Dxy Dxz \ 
Dxy Dyy Dy, (l) 
Dxz Dyz Dzz J 

where Dxx,Dyy,Dzz relate the diffusion flows to the concentration gradients in 
the X, y and z directions. The off-diagonal terms reflect the correlation between 
diffusion flows and concentration gradients in orthogonal directions. This diffusion 
tensor can be graphically represented as an ellipsoid. This ellipsoid takes the three 
eigenvectors of the matrix as principal axes (representing the three principal di- 
rections of diffusion). The length of the axes, related to the intensities of diffusion 
along them, is determined by the eigenvalues. Diffusion tensor images can thus be 
viewed as fields of ellipsoids. 

Classical image processing methods have been developed for scalar fields. As 
a result, early processing of DTI data first converted the tensor information into 
scalar data, f or instance focusing on th e scalar measure of fractional anisotropy 
(FA), see e.g. ( Alexander and GeelJ2QQ"Q) . However, the tensor nature of DTI data 



soon motivated a generalization of signal processing methodological frameworks 
to tensor fields. Methods based on the Riem annian geom etry of sy mmetric positive 
defin i te matrices have emerged in particular (IPennec et al, 2006: Fletcher and Joshii 
20071: ICastano-Moraga et aj l2006l : iGur and Sochenl . 120071: iBatchelor et aj l2005l : 



Lenglet et a i I2006l l2009ir because the geometric framework provides a nonlinear 



generalization of calculus in linear spaces. 

The present paper adopts the Riemannian framework to propose a novel metric 
and similarity measure between diffusion tensors. Our objective is to retain the 
geometry and computational efficiency of existing methods while addressing their 
main limitation, i.e. the degradation of the anisotropy information through aver- 
aging and interpolation. 

Our proposed similarity measure provides a concept of averaging and a concept 
of anisotropy that have the remarkable property to commute: the anisotropy of 
the average tensor is the average of the anisotropics. In that sense, our similarity 
measure is anisotropy-preserving. 

To arrive at this desirable property, we need a metric that compares separately 
the orientation of the tensors and their spectrum. Earlier methods based on this 
spectral decomposition suffer a strong computational obstacle that we remove 
thanks to quaternions. The use of quaternions to compare orientations is analog 
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to the use of Log-Euclidean computation to compare positive definite matrices. 
In that sense, our spectral quaternions framework for computation is analog to 
the Log-Euclidean f ramework for computation in the afiine invariant geometry 
(|Arsignv et ail2QQ7h . with similar computational gains. 

The paper is organized as follows: Section 2 reviews the geometry of diffusion 
tensors and introduces the novel metric. Section 3 develops the computational 
framework, that is a similarity measure that captures the desired geometry while 
being computationally tractable. Section 4 introduces the desired concepts of mean 
and anisotropy, which are shown to commute. Section 5 illustrates the properties 
of the proposed framework. Section 6 contains concluding remarks. 



2 Riemannian metric for diffusion tensor processing 

Even though the statistical analysis of diffusion tensor ima ges still largely relies 
on sc alar quantities such as the fractional anisotropy (FA) ([Alexander and Ged , 
l200Qh , the potential benefit of properly exploiting the tensor information of these 
data has been recognized early. The most common matrix representation of dif- 
fusion tensor is by elements of S^(3), the set of 3 x 3 symmetric positive definite 
matrices. 

Because this set is a nonlinear space, basic processing operation such as averag- 
ing and interpolation cannot be performed in the usual (euclidean) way. Defining 
and computing such quantities in nonlinear spaces equipped with a Riemannian 
geometry is a topic of active ongoing research, motivat ed by many applications 
in statistical signal processing (Smith, 2005; P ennec et al, 2006; Ando et al, 2004" 
Petz and Temesil . l2005HMoakher and ZeraJ. boilHMoakherl . l2005l : iBurbea and Raol . 
1982 : ISkovgaardi 119841 ). Focusing on S^(3), we briefly review the importance of 
the afhne-invariant metric and we introduce a novel metric. 



2.1 The afiine invariant Riemannian geometry of positive definite matrices 



A Rie mann ian framework for DTI pro cessing was first introduced in ( Pennec et aj 



I2OO6I) and ([Fletcher and Joshil l2007h . It is based on the parametrization of the 
space of symmetric definite positive matrices S^(3) as a quotient space, i.e. 

S+(3) ^ Gl(3)/0(3) (2) 

where Gl(3) is the space of general linear matrices (representing all the possible 
afiine transformations) while 0(3) is the space of orthogonal matrices. In matrix 
term, any element S G S^(3) can be represented by an invertible matrix A G Gl(3) 
through the factorization S = AA"^, but this representation is unique only up to 
multiplication by an orthogonal matrix U G 0(3) bec ause AA-^ = (AU)(AU)-^. 
The quotient geometry ([2]) is reductive ([SmitU [20051 ), which owes to the decom- 
position of an arbitrary 3x3 matrix (a tangent vector to Gl(3)) as the sum of a 
skew symmetric matrix (a tangent vector to 0(3)) and of a symmetric matrix (a 
tangent vector to S^(3)). As a consequence, the set S^(3) can be equipped with 
the special metric 

5s(e,r;)=tr(S-i/2^S-Ss-'/') (3) 
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which is called affine-invariant because of the property 

5gsg^ {G^G^, GriG^) = gsit n), VG e Gl(3) . 

The resulting Riemannian distance has the explicit expression 
d(Si,S2) = ||log(S-^/"S2Si-i/2||^ 

and inherits the affine invariance property of the metric, that is 

d(GSiG^, GSsG^) = (i(Si, S2) , VG G Gl(3) . (4) 

The metric invariance has two concrete consequences of importance for DTI pro- 
cessing. First, it is invariant to geometric transformation of the data such as ro- 
tation {i.e independent of the (arbitrary) orientation of the laboratory axes) or 
scaling (independent of the choice of length units), a very desirable property that 
respects the physical nature of MRI data. Second, the distance to identity 

d(S,I) = IllogSII (5) 

= /^log^A, (6) 



is indeed a measure of anisotropy, that is, does not depend of the orientation 
of the tensor S but only on its spectrum. The affine inva riant framewo r k has 
been used in several works in the context of DTI processing ( Fletcher et aj 2QQ9I; 



Castaho-Moraga et a ur and Sochenlj2007[|Batchelor et ail2Q05l:lLenglet et aj 

20Q6ll2QQ9[ ). 



2.2 A novel Riemannian metric based on the spectral decomposition 

A main motivation for the present paper is to address a limitation of the affine 
invariant metric graphically illustrated in Figure [TJ top: the anisotropy of ten- 
sors tends to be smoothened and reduced through averaging. If anisotropy is a 
key information of tensor, this means a loss of information that can quickly accu- 
mulate through the many averaging operations involved in data processing. This 
limitation, intrinsically related to the affine invariance of the metric ([3]), has mo- 
tivated alternative similarity meas ures that compare orientation and spectral in- 
formation separat ely ( Weldeselassi e et al , l20Q9t iTschumperle and Derichj , I2QQII : 
llngalhalikar et ail2QlQl ), starting from the spectral decomposition 

S = UylU^, (7) 

where U is an orthogonal matrix containing the eigenvectors of the tensor and A 
is a diagonal matrix whose entries are the eigenvalues of the tensor. Since tensors 
are symmetric positive definite matrices, A G D^(3), the space of diagonal matrix 
with positive elements, and the matrix U belongs to the orthogonal group 0(3). 
In this work, we will order the eigenvalues of the tensors such that Ai > A2 > A3 
and we will impose det(U) = 1, i.e., U belongs to the space of special orthogonal 
matrices, SO (3). 



An anisotropy preserving metric for DTI processing 



5 



\ • / 




Fig. 1 Mean of two cigar-shaped diffusion tensors. Top: Afiine-Invariant mean (jPennec et all , 
120061 ), bottom: mean computed with the method proposed in this paper, the spectral quater- 
nions method. 



Figure [TJ bottom, shows an example of a mean computed using the method 
developed in this work, and compares it to the affine-invariant mean (top). Both 
methods preserve the determinants of tensors (linked to the volume of ellipsoids) 
but, in contrast to the affine-invariant mean, the proposed mean conserves the 
shape of the tensor. This difference is of importance in the context of Diffusion 
Tensor Imaging, as the anisotropy property is a key information for instance in 
tractography. The following subsection will introduce in more details the proposed 
metric. 



A Riemannian metric based on the spectral d ecomposition Q is obtained as 



follows. Our construction builds upon earlier work (^Bonnabel and Sep ulchre! . [2QQ9I ) 
that aims at extending the affine invariant metric to 'fiat' ellipsoids, that is to pos- 
itive semidefinite matrices of given rank. 

Let S*^(3), be the subset of S^(3) where all the eigenvalues are distinct, i.e., 
< A3 < A2 < Ai. Each tensor S of S*^(3) can be identified to a diffusion 
ellipsoid whose equation is {x eM.^ S~-^x = 1}. Any such ellipsoid is determined 
by three orthogonal axes, and three positive values corresponding to the length 
of the ellipsoid's axes. Note that the axes are not oriented. Thus, the following 
identification holds 

S* + (3) ^:^T(3) X D+(3) 

where T(3) is defined the following way: each element of T(3) is a set of three or- 
thogonal directions in M"^, i.e., T(3) = {(si^i, S2t'2, S3V3) | Vi.vj = ^ij , (si, S2, S3) G 
MJ}. D^(3) is the set of 3 x 3 diagonal matrices with positive entries. When two 
or more eigenvalues are equal, the tensor is described by an infinity of elements 
of T(3). This case almost never occurs in practice and is thus excluded. The set 
T(3) can be advantageously parameterized by matrices of SO (3) based on the 
identification T(3) ^ S0(3)/G where G is the discrete group acting of SO (3) by 
multiplying two arbitrary columns jointly by —1 (i.e. G is the group of the three 
rotations of angle tt around the axes of the reference frame, plus the identity rota- 
tion). More prosaically, this identification expresses the fact that the orientation 
information contained in each tensor of S*^(3) can be described by 4 distinct 
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rotation matrices, see Figure [H As a result this set can be viewed as a product 
S* + (3) ^ (S0(3)/G) X D+(3). The space S0(3)/G is locally isomorphic to S0(3) 
(we say SO (3) is a quadruple cover of this set in the topological sense) and in- 
herits the usual metric gso(3) of ^^is latter space. We thus propose to define a 
Riemannian metric on S*^(3) expressed on the tangent space at 5* = UAU^ as 
follows 

fc(^)5SO(3) + 

where ^so(3) is the metric of the space S0(3) and 5'd+(3) ^^le metric of D+(3) to 
be defined later, and where /c is a weighting factor. 

The choice of the weighting factor k is motivated by the following fact: the more 




Fig. 2 Non-uniqueness of the spectral decomposition. Four different rotation matrices 
parametrize the same tensor. 



isotropic a tensor, the larger the uncertainty about its orientation (I Parker et aj 
^003). On the contrary, if two very anisotropic tensors are compared, the infor- 
mation contained in their orientation is very precise and the discrepency between 
orientations must be emphasized. We thus choose the factor /c as a function of 
the degree of anisotropy of the tensor. We let k be equal to zero when all the 
eigenvalues are equal, and tend to 1 as the anisotropy becomes large. 



3 Computationally tractable similarity measures 

Similarity measures for DTI should not only be geometrically meaningful. They 
must also be computationally tractable to scale with the combinatorial number of 
averaging operations involved in the statistical processing of 3D brain images. 

Computational tractability motivated a simplification of t he affine invariant 
metric known as the Log-Euclidean metric (jArsignv et all , [20071 ). Due to its numer- 
ical efficiency, the Log-Euclidean metric has been used extensively in recent stud- 
ies ("Coodlett et al", "2009; "Chiang et al", "2008"; "ingalhalikar et al", "2010"; "Castro et al', 
2007; Weldeselassie and Hamarneh, 2007; Arsignv et al, 2006; Fillard et al, 2007; 
Rodrigues et all . 120091: 1 Yeo et all . 120091 120081: iLepore et all . l2006l:lAwate et all . l2007l) . 
The geometry of the Log-Euclidean metric simply exploits the property that the 
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nonlinear space S+(3) is mapped to the linear space R^^^ by the (matrix) log 
mapping S logS. The Log-Euclidean distance is thus 

^iLE(Sl,S2) = ||log(Si) -log(S2)||F (8) 

where ||A||f is the Frobenius norm of matrices. This distance looses the affine 
invariance property ([4j) but retains the two concrete properties singled out in the 
previous section: (i) invariance to rotations and scaling and (ii) spectral nature 
of the distance to identity ([6]). It suffers the same averaging limitation as the 
affine-invariant metric illustrated in Figure [TJ 

In an analog way, we will exploit the euclidean embedding of rotation matrices 
in the linear space of quaternions to arrive at a computationally tractable similarity 
measure. 

The proposed Riemannian metric provides an infinitesimal notion of length 
on the set S*^(3), which allows to measure the length of any curve in S*^(3). 
The geodesic distance between two tensors is defined as the infimum of the length 
of curves joining those two tensors. Instead of calculating it, which can prove 
quite complicated and computationally expensive, we propose to directly define a 
similarity measure between tensors that approximates the geodesic distance, which 
is given in explicit form by 

diD(Si,S2) = ^(^1,^2) 4o(3)(Ul,U2)+4+(3) (^1,^2) (9) 

where dso(3)(Ui, U2) is a certain distance between two rotation matrices in the 
space S0(3), Ui and U2 are two well chosen representatives, while dD+(^^^{Ai, A2) 
is a certain distance between two positive diagonal matrices, both supposed to 
approximate the geodesic distances in the respective spaces SO (3) and D^(3). 
The function k is supposed to account for the anisotropics of Si and S2. 

Comparing intensities of diffusion 

As explained above, the similarity measure should be invariant under scalings and 
rotations. This implies that the distance between eigenvalues should be invariant 
to scalings (this distance is not affected by rotation). This is indeed the case as 
soon as their distance is of the form 

dD+(3)(^l,/l2) = ^^log2 (^^j (10) 

where Ai^^ and A2,i are the eigenvalues of Si and S2, respectively. 

With the choice (p^ , the distance ([9]) to identity coincides with ([6]), provided 
that the weighting factor k tends to zero whenever Si or S2 becomes isotropic. 

Comparing orientation of diffusion 

A longstanding computational problem of the spectral approach for processing 
of tensor images is the overparameterization of T(3) by orthogonal matrices (see 
Figure [2]), as we have seen that T(3) admits a q uadruple cov er by S^(3). In pre- 
ceding works using the spectral decomposition (jTschumperle and Derichd . I2QQ1I : 
IChefd' hotel et al , 120041 ), methods of realignment of rotation matrices have been 
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proposed to overcome this difficulty. Mathematically, if we denote by Ui the set of 
four rotation matrices representing the orientation of the tensor S^, these methods 
search for the pair of rotation matrices Ui and U2 which minimizes the distance 
between Ui and IA2 

d{UuU2) = ^min dso(3)(ULU2) (11) 

where Ui is the chosen reference in U\ and o?so(3)(Ui, U2) is the geodesic distance 
in the SO (3) manifold 

^SO(3)(Ul,U2) = ||log(ufU2)||. (12) 

This step can thus be expressed in our framework as the search for the geodesic 
distance in T(3) between two orientations in the sense of the natural metric of 
T(3), i.e., the metric inherited from its covering group SO (3). However, since this 
realignment is local and has to be performed for each tensor of any image, this 
procedure dramatically increases the computational cost of processing algorithms. 

It is thus the parametrization of diffusion tensors by rotation matrices that 
make the framework computationally expensive. But the orthogonal group is 
mapped to a linear space through the log operation, exactly as the space S^(3). 
This is the essence of the classical quaternion-based computational framework of 
rotations. 



Euclidean embedding of orientations 



The use of quaternions is very popular in robotics to efficiently compute over 
rotation matrices. The group of quaternions of norm 1, usually denoted Mi provides 
an universal cover of S+(3) and thus an universal cover of T(3). A unit quaternion 
is generally denoted by q = (a, V) where a is associated to the angle of rotation 
by ^ = 2arccos(a) and V is associated to the axis w of rotation through w = 
V/sin(^/2). From q, the associated rotation matrix R is given by 

(0 —wsO W2O \ 
W3O -wiO . (13) 
-W2O wiO / 

The construction of q from R is given by 

6 = arccos((trace(R) - 1) /2) (14) 



1 



2 sin 6' 



^3,2 — ^2,3 \ 

Ri,3-R3,i . (15) 



Finally, we have a = cos(^/2), V = sin(^/2)w. Note that the opposite quaternion 
given by (—a, —V) represents the same rotation matrix. Using this representation, 
rotations can be manipulated as Euclidean vectors, which decreases the compu- 
tational cost. For example, we can use a simple Euclidean norm to compute the 
distance between two unit quaternions 



G^(qi,q2) = ||qi - q2 



(16) 
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This is the so-called chordal distance on the sphere of quaternions, which is 
a good approximation of the geodesic distance of SO (3) for close enough quater- 
nions in Ml. This latter group being an octuple cover of T(3), processing tensor 
orientations with quaternions could seem more involved, but proves to be much 
more economical, computationally speaking. Indeed, let Qi be the set of eight 
quaternions representing the orientation of the tensor S^. We propose the follow- 
ing distance on T(3) 

d{Qi,Q2) = min ||q[ -q2|| (17) 

q2GQ2 

where q^ is one element of Qi chosen as a reference. This element can be chosen 
arbitrarily, and the distance is a good approximation of the geodesic distance on 
T(3), inherited by the geodesic distance on SO (3). Since qi and q2 are of unit 
norm, this distance can be advantageously expressed as 

||qi -qsiP =2-2qi.q2 (18) 

Using Eq.([T8]), the distance between Ui and U2 can be computed through 

q2 = max q[.q2 (19) 

q2GQ2 

d(Qi,Q2) = ||ql-q§|| (20) 

where q2 is called the realigned quaternion. 

Compared to (pJ^ , the computation of (p^ and (|2Q)) is very fast. Indeed, the 
eight scalar products qi.q2 can be computed through a single matrix product be- 
tween the 1x4 vector representing (qi)"^ and the 4x8 matrix formed by the eight 
quaternions q2. In contrast, computing the distance ()lip requires four logarithms 
of product of 3 X 3 matrices, which is expensive. The selection of the parametriza- 
tion of rotations as quaternions thus enables the framework to be computationally 
tractable. Moreover, the distance corresponds to a natural distance on T(3), and 
is thus invariant to rotations. 

The similarity measure ([9]) between two tensors Si and S2 becomes 

dio (Si , S2) = fc I |q^ - q§| |2 + ^ log2 f^) (21) 

where the function k of the anisotropics of both tensors will be exactly defined in 
the sequel (section [4]) . This measure will be hereafter called the 'spectral quater- 
nions' measure. 

Table [T] illustrates the numerical savings obtained by proper euclidean em- 
bedding of the desired Riemannian geometry. The saving of using quaternions as 
opposed to distances over rotation matrices is of the same order as the saving of 
using the Log-Euclidean metric instead of the affine-invariant metric. The table 
suggests that the numerical cost of the proposed similarity measure is competitive 
with the numerical cost of Log-Euclidean distance. 

Figure [3] shows the variation of the spectral quaternions distance and the Log- 
Euclidean one when the characteristics of the tensors are smoothly varied. In each 
subfigure, the distance between each tensor of the set and the 'central' tensor of 
the same set is computed. When only the eigenvalues are varied (left of the figure). 
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Table 1 Computational time of computing 1000 distances between a reference and random 
samples from a Wishart distribution. The computations are performed on a Intel Core 2 Duo 
2,66 GHz with 4Go of RAM machine using a (non optimized) MATLAB code. 



Affine- Log-Euclidean 


Spectral 


Spectral- 


invariant 




quaternions 


0,47 s 0,17 s 


0,65 s 


0,11 s 



both the Log-Euclidean and the spectral measure give the same linear variation. 
In the center subfigure, only the angle of the principal eigenvector is varied. In 
this case, the spectral measure varies smoothly, in sharp contrast with the Log- 
Euclidean metric, which is very sensitive to small angle variation. The two effects 
are combined in the right subfigure, when both eigenvalues and orientation are 
varied. 




Fig. 3 Evolution of the different similarity measures when the properties of the tensors are 
smoothly varied. Left Variation of eigenvalues. In this case, the spectral and Log-Euclidean 
measures give the same results. The distances vary linearly. Center Variation of the angle of the 
dominant eigenvector. While the variation of the spectral distance is smooth, this is not the case 
for the Log-Euclidean metric. This metric appears to be more sensitive to variations in rotation, 
while the spectral distance is robust to these variations. Right Variation of both eigenvalues 
and angle. The Log-Euclidean measure is not regular, while the spectral quaternions distance 
is almost linear. 



4 Anisotropy commutes with averaging 

In this section, we build upon the proposed similarity measure to define an anisotropy 
index and an averaging operation among tensors that have the remarkable prop- 
erty to commute: the anisotropy of the average is the average of the anisotropics. 
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4.1 Hilbert anisotropy 

Geometrically, any anisotropy s calar mea s ure sh ould be a scale invariant distance 
to identity. The Hilbert metric ( Koufanvi |2006[ ) is precisely a projective distance 



that can be defined in arbitrary cones. It leads to the following definition that we 
refer to as Hilbert anisotropy (HA) in the sequel: 

HA = dH(S,I)=log(^) (22) 

where Amax is the maximum eigenvalue of S and Amin, the minimum one. The HA 
index possesses all the required properties for an anisotropy index, i.e. 

— HA > and HA = only for isotropic tensors. 

- HA is invariant to rotations: HA(S) = HA(USU^) for all U G 0(3). 

- HA is invariant by scaling, HA(S) = HA(aS), \/a G IR+ (it means that anisotropy 
only depends on the shape of the tensor and not on its size). 

— HA is a dimensionless number. This property is desirable and natural, as the 
anisotropy of the tensor physically refiects the microscopic anisotropy of the 
tissues, which is independent from the diffusivity. 

Figure |4] illustrates a comparison of HA with three p opular anisotropy in- 
dices : fractional anisotropy (FA), rel ative anisotropy (RA) (iBasser and Pierpaolii 
[1996), the geodesic anisotropy fG A) ([Fletcher and Joshil l2QQ7h . This figure is con- 
structed with eigenval ues equal to t, (1 — t)/2, (1 — t)/2, t g]0, 1], as explained in 
([Batchelor et all , |2QQ5[ ). The diffusion tensor varies from planar to spherical for 
t g]0, 1/3] and then tends to become unidirectional when t increases. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Fig. 4 Evolution of four indices of anisotropy with respect to modifications of eigenvalues. The 
difference between 'linear' indices of anisotropy (FA, and RA, ^) and 'geometric' ones (HA, 
O and GA,D) i s clea rly visible. The eigenvalues are given by t, (1 — t)/2, (1 — t)/2, t G]0, 1], 
iBatchelor etal (I2005I V 



12 



A. Collard et al. 



We use the HA index to define a weighting factor k in (|2T]) . The empirical 
choice adopted in this paper and that can be customized by the user is 

MHA„HA.) = l+tanh(3HA,HA.-7) ^^3^ 

which smoothly varies from zero when one of the tensor is isotropic to one when 
the two tensors are strongly anisotropic. 



4.2 Weighted mean of Diffusion Tensors 



The Riemannian mean of a set of tensors is defined as the tensor that minimizes 
the total squared distance to all tensors (see ([Fletcher and Joshii l2Q07l )V Here we 
use our similarity measure to provide a numerically tractable approximation of 
this quantity. 

— Weighted mean of two tensors Consider two tensors Si and S2 with respective 
weights wi >0 and W2 > {wi -\- W2 = 1). 

- The eigenvalues of the mean tensor are defined as 

X^^i = exp(i(;i log Xi^i + W2 log A2,i) (24) 

which corresponds to a Riemannian mean on D^(3) for the logarithmic 
distance {i.e. a geometric mean of positive numbers). 

— The mean orientation is defined by means of the mean quaternion 



which is a weighted chordal mean in the space of quaternions. Note that 
qi and q2 in ([25]) must be realigned quaternions as in ([2Q|) . 
Weighted mean of more than two tensors Consider N tensors Si, S2, . . . , S^r 
{N > 2) with respective weights wi,W2, • • • w^, J2i '^i = 1- Using our similarity 
measure, we propose a method that provides a well-defined and commutative 
mean. 

— Similarly to the case of two tensors, the eigenvalues of the mean are defined 
as the geometric mean of all eigenvalues, i.e. 

N 

yl^ = exp(^^,log(A)) (26) 

i=l 

- To make the definition of the mean unique, we must select a quaternion, q^, 
which will serve as a reference for the realignment of all other quaternions. 
To this end, we first compute the weighting factors ki = /c(HA^, HA^) , Vi = 
1, . . . , A^. We then choose as a reference the tensor Sr for which the product 
Wrkr is the highest {i.e. the most informative tensor). The realignment pro- 
cedure is then applied to all the other tensors, and the realigned quaternions 
are denoted q^^^. 
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— The mean orientation is defined as 

ZJi 

The mean is commutative, because the reference quaternion does not 
depend on the labeHng order of the tensors. 

The chordal mean (|25]) , (|27|) of quaternions is the solution ( Dai et al |201Q[ ) of 



N 



R = arg min ^ ^.^t (R, Ri ) (28) 

i=l 

where dquat is the quaternion distance between rotations and are the rotations 
to be averaged. It is thus a Riemannian mean over the set of rotations for the 
chordal distance defined on the set of quaternions. Since the quaternion distance 
is a good approximation of the geodesic distance between rotations for small ro- 
tations, the quaternion mean is viewed as a good approximation of the real mean 
rotation. 



4.3 Properties of the mean 

The proposed spectral-quaternion mean (denoted in the following) enjoys two 
'information- preserving' properties of interest for DTI processing. 



Determinant ofS^ From Eq. ()24|) . one easily shows that 

det(S^) = ey:p{wi log(det(Si)) + W2 log(det(S2))) 

i.e. the determinant is the geometric mean of the determinants of the tensors. 
The same property holds in the case of > 2 tensors, thanks to Eq. ([26]) . 



N 

det{S^) = exp(^i/;i log(det(Si))) 

i=l 

The Log-Euclidean mean shares the same property (jArsignv et all . l20Q7h . This 
ensures that there is no swelling effect (increase of the volume of the ellipsoid) 
during the averaging of tensors. 

Anisotropy of The anisotropy of the mean is the arithmetic mean of the 
anisotropics of the tensors. In other words, anisotropy commutes with averaging 
of tensors. Indeed, if HAi is the anisotropy of Si and HA2 the one of S2, we 
have that 



HA^ = log 



wi log ( ) + '^2 log ^^'^ 



.Al,3 / \A2,3 

= wiYLAi + i6;2HA2 
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In the case of more than two tensors, it can easily be shown that 

N 

HA^ = 

i=l 

This important property is not satisfied by the Log-Eucfidean method, neither 
by the affine-invariant mean. It is because of this property that the proposed 
metric addresses the shortcomings of the affine invariant metric pointed in 
Figure [TJ 



5 Information preserving interpolation 

As a direct application of the weighted means introduced in the preceding section, 
we propose an anisotropy preserving interpolation scheme for diffusion tensor im- 
ages. Other examples of application of weighted means include the computation 
of statistics about group of tensors. 



As already mentioned in ( Zhang et allJ2QQ6l : lKindlmann et allJ20Q7l : lArsignv et all . 



l2QQ7l) . an adequate interpolation method is important for the processing of diffu- 
sion tensor images and particularly for the extension of usual registration tech- 
niques (for scalar images) to the case of tensor images. This interpolation scheme is 
necessary to resample images. Here, we provide a direct generalization of classical 
interpolation method, where the interpolated value is computed as the weighted 
mean of the original tensors. Equations ([24|) and (|25)) directly provide interpolating 
curves in the space of diffusion tensors, with, if t is the parameter of the inter- 
polation, wi = (1 — t) and W2 = t. Those curves are reasonab le approximation of 
geodesic curves employed in Riemannian interpolation theory (jPletcher and Joshil , 
l2Q07l ). 

The scheme for the bi and tri-linear interpolations of scalars can not been ex- 
tended to tensors, because unlike linear interpolations, interpolation along curves 
does not commute. This is why the commonly used solutio n is to compute in- 



terpolation through a weighted average of d iffusion tensors ( Pennec et aj l2QQ6l : 
IPletcher and Joshi l2Q07l : lArsignv et ail2QQ6l) . The weight associated to each ten- 
sor is a function of the grid distance between this tensor and the location of the 
interpolated tensor. In this work, if (xi,X2,X3) G [0,1] x [0,1] x [0,1] are the co- 
ordinates of the interpolated tensor and (ai, 0^2,0^3) ^ {0, 1} x {0, 1} x {0, 1} the 
coordinates of the point a of the grid, the following function will be used 

3 

Wa{xi,X2,X3) = ]^(1 - tti + {-iy~^'Xi). 
i=l 

Pigure [5] shows the curve interpolation between two tens ors using bot h the Log- 
Euclidean and the spectral quaternions frameworks. As in (|Zhoul , [2"oiO h, the varia- 
tion of the main information conveyed by the tensors is also shown. As previously 
shown, the Hilbert anisotropy is linearly interpolated by the novel framework, 
while this information is significantly degraded in the Log-Euclidean method. A 
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similar behavior is found for the evolution of the fractional anisotropy. Both meth- 
ods geometrically interpolate the determinant. It is also interesting to analyse the 
difference in 0, the angle between the first eigenvector of the first tensor and the 
first eigenvector of the weighted mean. While the spectral measure produces a 
quasi linear interpolation of this angle, this is not the case for the Log-Euclidean 
framework. 



• • • 9 i I I I I 



* ^ ^ f f ( ( I 




Fig. 5 Geodesic interpolation between two tensors. Top: Log-Euclidean interpolation. Bottom: 
Spectral-quaternion interpolation. Clear differences can be seen between these two methods, 
namely for the evolutions of degrees of anisotropy. 



Eigure [6] illustrates the invariance by rotation of the curve interpolation of the 
spectral quaternions framework. Although the interpolations are different (in par- 
ticular for the shape of the tensors), the Log-Euclidean framework also possesses 
this property. 

Using the method described above for computing the weighted means of many 
tensors, the interpolation of four tensors at the corners of a grid can be computed, 
as illustrated in Eigure [71 where colors of the tensors is determined by HA. The 
four tensors at the corners have equal Hilbert Anisotropy. As a consequence, all 
tensors interpolated by the spectral method conserve the same anisotropy, while 
this is not the case with the Log-Euclidean method. 



6 Conclusions 

In this paper, we have introduced a novel geometric framework for the processing 
of Diffusion Tensor Images. This framework is based on the spectral decomposition 
of tensors. The main advantage of this method is to preserve the anisotropy of ten- 
sors during the processing. Moreover, it possesses all the important properties of 
existing metrics, such as the invariances and the preservation of other information 
as the determinant and the orientation of tensors. The preservation of anisotropy 
is rooted in the spectral decomposition of the tensors, which allows for comparing 
separately rotations and the spectral information. 
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Fig. 6 Geodesic interpolation between two tensors, illustration of the invariance by rotation. 
Left: Log-Euclidean interpolation. Right: Spectral quaternions interpolation. The colors of 
tensors are based on the anisotropy. From top to bottom, the original tensors are rotated. 
With both frameworks, it is clear that the interpolated tensors are rotated in the same way, 
which means these frameworks are invariant by rotation. However, a clear difference can be 
seen between the frameworks, since the spectral quaternions one preserves the anisotropy, while 
the Log-Euclidean one produces a degradation of this information. 



n 



Fig. 7 Geodesic interpolation between four tensors at the corners of the grids. Lejt: Log- 
Euclidean interpolation. Right: Spectral interpolation. Colors of the ellipsoids indicates their 
anisotropy from red (high anisotropy) to yellow for smaller anisotropics. While the anisotropy 
is equal for each tensor computed by the spectral quaternions method, this is not the case with 
the Log-Euclidean one, where a decrease of the anisotropy can be seen. 



Computational obstacle s previ ously faced in similar attempts ( Tschumperle and Derich3 , 
I2OOII: I Chefd'hotel et all , |200J) are circumvented by embedding the set of rota- 
tion matrices in the space of quaternions, long used for its numerical efficiency in 
robotics. 

The resulting interpolation method retains the computational tractability and the 
geometry of the Log-Euclidean framework but addresses a limitation of this frame- 
work regarding the degradation of anisotropy. 

Although several illustrations of the paper exemplify the potential benefit of pre- 
serving anisotropy through averaging and interpolation operation encountered in 
statistical process, the benefits of the proposed framework remain to be demon- 
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strated on real data. Registration and tractography are two particular areas where 
the advantages of the proposed method should be evaluated quantitatively. 
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